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Abstract 



We present a preconditioning method for the multi-dimensional Helmholtz equation 
with smoothly varying coefficient. The method is based on a frame of functions, that 
approximately separates components associated with different singular values of the op- 
erator. For the small singular values, corresponding to propagating waves, the frame 
functions are constructed using ray-theory. A series of 2-D numerical experiments demon- 
strates that the number of iterations required for convergence is small and independent 
of the frequency In this sense the method is optimal. 
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1 Introduction 

In this paper we will describe a preconditioning method for a class of Helmholtz operators 



where c{x) and a(x) are positive, smoothly varying coefhcents, d is the frequency parameter, 
and L is a typical length of the domain, which is inserted to make a dimensionless. See 
dU HH [231 EH 021 02] for examples of preconditioners from the literature. 

For elliptic equations there are good preconditioning methods, such as multigrid and 
wavelet diagonal preconditioning [6j [8] . Their effect can be explained using analysis in the 
position-wave number domain, the so called phase space. With this point of view we first give 
a short discussion of wavelet diagional preconditioning, and then turn to its generalization to 
the Helmholtz equation. 

Consider, as an example problem an elliptic partial differential equation Bu = /, where 
B is an invertible, symmetric, second order elliptic operator on a perioc domain (we will not 
go into the issues related to boundary conditions.) The main ingredients of the method are 
a basis transformation F, given by a wavelet transform, and an invertible diagonal scaling 





matrix W. The matrix F* is a wavelet reconstruction operator, its columns are the elements 
of the wavelet basis, denoted by F is a wavelet decomposition operator. The elements of 
the diagonal of W satisfy ~ 2 _2 l /i l, where \fx\ denotes the scale index of the wavelet, and 
by convention = — 1 corresponds to the scaling functions. A symmetrically preconditioned 
operator B is then example given by an expression of the form 

W 1/2 FBF*W 1/2 . (2) 

Left- and right preconditioned variants also exist, they are given by WFB and BF*W . 

The triple (F*,W~ l ,F*) can be viewed as an approximate SVD. A true SVD (U, S, V) 
would yield a perfect preconditioner, given by S^^U'BVS- 1 / 2 = 1. Equation (B yields 
an approximation to the identity in the sense that W 1 / 2 FBF*W 1 / 2 is boundedly invertible, 
uniformly in parameters such as the grid constant. 

The singular values of partial differential operators (PDO's) with constant coefficients 
follow from Fourier analysis. The symbol of the operator is defined by 

= e - ix <(Be ix <). (3) 

The singular functions are complex plane waves and the singular values -B(£ M ), with the 
wave vector. Wavelets are localized in the Fourier domain around ||£|| = 7r2' M ', away from 
zero and infinity, according to the vanishing moments and smoothness properties. Here being 
localized means that its Fourier transform decays rapidly away from some bounded region. 
Using this it can be shown that they act indeed as approximate singular function with value 
~7r 2 2 2 W 0[9]. 

For variable coefficient PDO's, or pseudodifferential operators, the analysis becomes more 
complicated, but in many cases similar results hold. The symbol becomes a function of (x, £) 
[111211 151]. If the symbol doesn't vary too much over the area in phase space where a function 
is large, then one can often still show that the operator acts as an approximate multiplication 
by the value of its symbol at the phase space support, or is an approximate singular function. 
This idea is well known in semiclassical analysis, see e.g. [281 [25l [T71 118j . and also occurs in 
numerical analysis |32j . 

With a wavelet basis is associated a tiling of phase space. The phase space can be divided 
into disjoint sets 7^ that approximately correspond to the support of the wavelets. In 1-D, the 
multi-index // has two components, (j, k), and the tile is given by 2~ 3 [k, k + 1] x 2 J ([ir, 2ir] U 
[— 2-7T, 7r]) for j > and by [k,k + 1] x [— ir, ir] for j = —1. The tiling and the matrix W are 
adapted to the operator B in the following sense 

C- 1 < w^a B {x,0 < C, for all € 7^, for all /x. (4) 

In words: The factor u;^ compensates for the action of the operator, which is an approximate 
multiplication, assuming the symbol doesn't vary too much over the tile. This explains the 
choice of the wavelets and the matrix W, even if it cannot be taken literally and must be 
supplemented with appropriate technical arguments. The property Q will serve as the guiding 
criterion for the design of a preconditioner for the Helmholtz operator. 

Applying the criterion Q to ([!]) , one finds that the requirement that F is a basis transform 
is very restrictive. To alleviate this requirement, methods based on frames have been proposed 
|30| [21 [20l I19j. If F is a tight frame with F*F = I, one can use a symmetric preconditioning 
of the form 

F*W 1/2 FBF*W 1/2 F (5) 
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Left- and right-preconditioned variants are given by 

BF*WF, and F*WFB. (6) 

A second form of left-preconditioning is given by WFB. The latter form leads to an overde- 
termined system of equations that is to be solved in the least-squares sense. The operators 
(§ and (§ are invertible if the diagonal entries of W are real and strictly positive. 

The criterion Q has important consequences for the preconditioning of the Helmholtz 
equation. Let's start with the situation that the coefficients c and a are constant. Clearly 
a wavelet tiling can't be used in a neighborhood of the set ||£|| ~ — . The wavelet will 
be of size 0(u>) in the radial direction, and cover a wide range of values of the symbol from 
= 0(u) to |i?(£)| = 0(oj 2 ). Based on the criterion they can only be used some distance 
0(u}) away from this set. Therefore, the generalization of wavelet diagonal preconditioning to 
the Helmholtz equation must use a different set of basis or frame functions, with a different 
phase space tiling. 

The tiling must be such that the tiles around ||£|| = ^ must be of size at most O(l) in 
the radial Fourier direction. It follows that they are of size 0(1) in the corresponding spatial 
direction. In other words, the frame functions have a macroscopic spatial extent, also when 
uj increases and the wavelength becomes smaller. This conclusion can also be reached from 
a completely different way, by considering the flow of information in an iterative algorithm. 
The distance the information can travel is proportional to the size of the support of the 
frame functions. Since we aim at convergence in a number of steps independent of u, and 
since the solutions of Hu = f contain waves propagating over 0(1) distance, there must be 
frame functions with O(l) size of the support, and it is natural that these correspond to the 
propagating wavenumbers with ||£|| ~ — . One can also observe that the number of frame 
functions with small approximate singular values is much larger than for the elliptic case, it is 
on the order of O(l) time the area of the set ||£|| = ^ . If the number of points per wavelength is 
kept constant while u may vary, this leads to 0(N d ~ 1 ) basisfunctions with small eigenvalues, 
where N is the number of grid points in each direction, and d is the dimension. For constant 
coefficients of course a Fourier basis can be used. 

In case of variable coefficients the situation becomes more complicated. The tiles with 
small value of the symbol \H(x, £)| must satisfy ||£|| ~ over an O(l) domain. The 
wavenumber content of the basisfunction must depend on the position. The frame functions 
must therefore be adaptively chosen, depending on c and uj. This is unlike many existing phase 
space transforms |24| . and also unlike certain transforms that have been used for hyperbolic 
PDE (but not for preconditioning) ESI E] . 

In this paper we show that it is possible to construct such a frame. The construction uses 
the WKB approximation in combination with a paraxial approximation for the wave equation 
(because ray theory is used it can be compared with the work of Brandt and Livshits [SJ [23] 
except that they treated the variable coefficient case only in one dimension.) Our second main 
result concerns the convergence of the preconditioned operator. Using an implementation in 
Matlab we find that the number of steps required to converge becomes bounded by a small 
number, independent of the frequency. 

While the number of steps is small, the cost per step is still relatively high. Further 
research into these transforms might improve this, as could the fast algorithm of [3] for the 
3-D case. We find that larger problems can be solved than with the direct method because 
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the memory requirement is reduced. The computation is split in two steps, a preparation step 
and an execution step, such that the results of the preparation step can be reused for each 
right hand side. We find that the computation time for the execution step is lightly increased 
compared to the direct method. The smoothness requirement on the medium, the inclusion 
of other boundary conditions or a damping layer, and the extension to the 3-D case are topics 
for further research. 

The remainder of the paper is structured as follows. In section [2] some properties of our 
class of Helmholtz operators are discussed. The construction of the preconditioner in the 
continuous setting is the topic of sections [3] and |4} Some aspects of the discretization and 
implementation are then explained in section [5] Section [6] shows the results of our numerical 
experiments. We end with a short discussion. 



2 A class of Helmholtz operators 

The class of Helmholtz operator we study is given in 

0. We assume that there are constants 

C\ and C2 such that 

< Ci < c(x) < C 2 . (7) 

The imaginary part of H is responsible for damping. We assume a{x) is O(l) throughout the 
domain. So there must be C > such that 

i < a(x) < C. (8) 

We assume the domain is rectangular. For simplicity, the boundary conditions are assumed 
to be periodic. The spatially varying wave length is given by 2n< ^ x ) , The situation of interest 
is when the domain is many wave lengths large, but not so large that a discretization of the 
Helmholtz equation becomes impossible. 

With the periodic boundary conditions the part of H without the imaginary coefficient, 
which we denote by 

is real and selfadjoint. It has real eigenvalues that may be close to zero. The imaginary part 
influences the singular values of the operators H and A. It leads to the lower bound 

\\ Hu \\ > J » ( 10 ) 

where a m i n is the minimum of a over the domain, and c max the maximum of c. 

We consider a finite difference approximation A of H based on the standard five point 
stencil for the Laplacian 

1 ( to 2 icti jU)\ 

{Au)ij = (4uij - u i+ ij - Ui-ij - Uij + i - Uij-i) + I ~^2~ + Lc J u ^r ( n ) 

The preconditioning method is based on the properties of the continuous operator, therefore 
we expect the results to be valid for other discretizations as well, as long as they approximate 
the continuous problem with reasonable accuracy, see also the remarks later on. Setting on j 
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to zero in (11) yields the finite difference approximation A\ of H\. The matrix A\ is real and 
selfadjoint, like Hi, so that A satisfies 



Pu|| > (12) 

Based on the largest eigenvalue of the discrete Laplacian, given by 8h~ 2 , we find the following 
estimate for A 

8 uj 2 uj 
+ i- 



\\Au\\ < 



|u||. (13) 



The number of points in the domain, and the grid size h are in practice often chosen based 
on the number of wavelengths that fit in the domain. For example using a rule to keep a 
certain number of points, say N w , per wavelength. Reasonable values are of the order of 15 
or 20. This means 

h= 2n ^. (14) 

iVwW 



dm ax 



Equations (12) and (13), and the assumption that a is not too large, so that 

< A imply the following bound for the condition number 

K{A) < 2 f^ C 2 m " ( 15 ) 

n a minC m i n 

In other words IC(A) ~ uj if uj — > oo with the given choice of h. 

We have made the particular choice that the damping term has coefficient 0(uj). This 
leads to the suppression of certain resonance effects, where small eigenvalues occurring around 
specific frequencies lead to numerical difficulties. At the same time, the solution to Helmholtz 
equation still consist of propagating waves with propagation distance on the order of the 
domain length. A damping term 0{uj 2 ) would lead to a propagation distance on the order of a 
fixed number of wavelengths when uj — > oo. The problem would become essentially elliptic. For 
elliptic equations, wavelet or multigrid preconditioning methods lead to a uniformly bounded 
operator, i.e. the problem of finding a preconditioner is essentially solved. The resonance 
effect that were just mentioned depend on subtle global properties of the medium. Some 
of the small singular values can likely be removed by using outgoing radiation boundary 
conditions, but it appears that internal wave patterns can also lead to large values for the 
norm of the resolvent. In this respect it could be interesting to study consequences of the 
results from the semiclassical analysis of the Schrodinger operator and its resonances, see e.g. 
[16 1 . It is clear that the ideas of this paper do not take into account the global properties 
referred to, so we make no claims concerning the situations without damping present. 



3 Frame transform adapted to the Helmholtz operator in one 
dimension 

In this section we design a frame transform adapted to the Helmholtz operator with the domain 
[0, 1] and periodic boundary conditions. The frame transform will consist of two steps. First 
we separate the function into three components that are essentially supported in three regions 
of phase space: The large wave numbers, the positive wave numbers of wavelength scale, and 
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the negative wave numbers of wave length scale. The small wave numbers (smaller than 
wavelength scale) are treated together with the wavelength scale wave numbers. For each 
of these components a transform corresponding to a tiling is applied: Fourier tiling for the 
large wave numbers, and a modulated Fourier tiling (associated with a modulated Fourier 
transform) for the wave length scale wave numbers. The modulation parameter is a spatial 
function depening on the function c, and will cause the required curved space tiling. The idea 
of a modulated Fourier transform for preconditioning is new to our knowledge. 
We first apply a windowing operator 



(16) 



using 

The function Xj £ C 00 ^) to minimize to keep signals localized in space. They must satisfy 
Y^j=iXj = 1; then this map satisfies G*G = I, so that a tight frame can be constructed. 
The function xi wm be supported away from the zero set of the real part of the symbol, 
we require xi(0 = on [— 1.5w/c m i n , 1.5w/c m i n ]. The function X2 will be supported around 
the positive zero set of H\. We require X2(£) = outside [— 0.4w/c max , 2w/c m i n ]. Similarly, 
the function %3 will be supported around the negative zero set of H\ , and X3 (£) = outside 
[-2w/c m in,0.4o;/c max ]. Within the overlap regions [-2w/c min , -1.5w/c min ], [-0.4w/c max , 0.4w/c n 
[1.5w/c m i n , 2w/c m i n ] we use a dilated and translated version of a smooth cutoff function h(x), 
such that h(l — x) 2 + h(x) 2 = 1, and given by [2] 



h(x) = < 



( if x < 0, 

if < x < 1, (18) 



2 — & 
e x +e 



1 if x > 1. 

For each of the Uj we use a basis transform. For u\ we use the standard Fourier transform 

Vj = <'- ri - (19) 
The operator —d 2 — uj 2 /c(x) 2 acts on this as 

P lVj = _ oo 2 /c{x) 2 ) Vj (20) 

The preconditioning factor is 



(eJ-^VcLan) 2 +(;^n ' (21) 

For U2 we use the functions 

Vj = e ^T(x)+ixi 3 ^ 

with T(x) = Jq c(s) _1 ds (the travel time), or with a small linear term added, so that w(T(l) — 
T(0))/(27r) G Z and e* wT hence satisfies the periodic boundary conditions. The values of £j 
are £j = 2irj. The addition of the term iuT(x) causes an additive shift in the so called 
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Figure 1: Example of curved phase space tilings adapted to H for £ > 0. (a) Modulated 
Fourier tiling; (b) modulated wavelet tiling. 

instantaneous phase [23] by 

^l)x' rectangular tiles associated with the Fourier modes 

v - e lx ii hence receive an x-dependent shift in the vertical direction of the (x, £) plane. The 
tiling associated with this basis is given by 

0<x<l, " +£ +7r j €Z . (23) 

c(x) c(x) 

We define the scaling factor for preconditioning by 

because we find that the value of the symbol can be uniformly estimated by w^ 1 on the 
intersection of the support of v,2 and the domain of the tile. The tiling for an example c(x) is 
displayed in Figure [Tj We treat K3 similarly with 

Vj = e-^W+^i, (25) 



and preconditioning factor given again by (24). Denoting by F the full transform, it follows 



that F*F = I, i.e. the transform corresponds to a tight frame. 



Instead of (22) and (25) it is also possible to use modulated wavelets, see the tiling in Fig- 
ure [TJ This can be attractive because a wavelet transform has cost O(N), versus 0(N log N) 
for the Fourier transform. 

Some computations using the above methods that were done early on in this research gave 
encouraging results. As the focus of this paper is on the multi-dimensional case we will not 
further report on this here. 

4 Frame transform adapted to the higher-dimensional Helmholtz 
operator 

In this section we devise a frame of functions, assocatiated with a tiling of phase space adapted 
to \H(x,£)\, in dimension d = 2 or higher. The construction is done in the continuous setting. 
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Discretization and implementation aspects are discussed in the next section. 

Like in the one-dimensional case the most challenging region of phase space is around the 
zero-set (characteristic set) of the real part Hi, given by 

Char^) = S.(x,0 € T d x R d ; ||£|| = ^y} • (26) 

Here T d denotes the (i-dimensional torus, also written as [0, l] d with periodic boundary con- 
ditions. In the Fourier space, the direction normal to Char(ffi), a tile has size at most 0(1), 
implying that the size in the corresponding spatial direction is at least 0(1). The tile must 
follow the curved surface of Char(f/i), a property not present in existing multi-dimensional 
transforms, see e.g. Mallat [21]. 

Just like the previous section we will not construct a single, global tiling. Instead the first 
step of the frame transform is a division into different regions of phase space. Regions asso- 
ciated with sufficiently large and sufficiently small wave numbers can be treated with Fourier 
preconditioning. The other regions are chosen such that the part of characteristic set they 
contain is a graph, say £j = f(x, £i, . . . , £j'-i> • • • , Cn)- We call Xj the preferred coordi- 
nate. We summarize the remainder of the construction. Locally in phase space the Helmholtz 
equation becomes equivalent to a first order evolution equation in the preferential direction. 
Using theory for hyperbolic pseudodifferential equations an abstract frame is constructed. To 
obtain a computable set of functions a WKB approximation is made, which finally yields the 
result. 

So the first step in the frame transform is the separation of u into components that live 
in different regions of phase space. Localization is done in the Fourier space in radius (scale) 
and angle, and in the position space to certain rotated rectangular regions. There are three 
wave numbers scales, the small wave numbers ||£|| < w/c max , the wave length scale wave 
numbers w/c max < ||£ll — w / c min, and the large wave numbers ||£|| > uj/c m \ n . We will use a 

~h\ 

map F sca ie :«(->■ lu 2 I using 

with three smooth cutoff functions Xj(£)> 3 = 1) 2, 3 to separate the small, wavelength scale and 
large wave numbers. Again the cutoff functions will satisfy ^ x| = 1 to have F* caie F scai \ e = I. 
By ki, &2 we denote constants to separate the low and wave length scale wave numbers. For 
||£|| going from k\ to k 2 Xi goes smoothly from one to zero and X2 smoothly from to 1 
according to the standard cutoff function h. By A;3,A;4 we denote constants to separate the 
wavelength scale and high wavenumbers. 

For further localization in the Fourier space, we will distinguish several directions in the 
Fourier space, each direction described by a unit vector a G S"^ -1 . Using cutoffs, the support 
in the Fourier domain will be restricted to angular wedges around each angle a. This angle 
will also determine the spatial localization. We consider new coordinates (y,z), with y in a 
direction and z = (zi, . . . , z^-i) normal to this direction. We cover the spatial domain with a 
small number of coordinate boxes in these coordinates. The new coordinates will have y = 
in the center of each box. Because of the periodic boundary conditions, in our 2-D numerical 
examples we chose periodic bands. To obtain this, a is set such that tan(ct) is a ratio of two 
small integers. 
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After phase space localization and coordinate transform, the preconditioning can be done 
according to the following symbol, for which we locally have P ~ H\ 



U) I uj 2 



P 2 Cmean \ V c(y,z) 2 ^ j ' 



Moreover, this can be further modified outside the region of phase space that the analysis was 
restricted to. Here this is useful to deal with the singularity of the square root at ||£ 

We define a C 2 function by 



S f 8z \ _ J V 1 - s l if M < sin O) (29) 

cos(a) — (\s z \ — sin(a)) tan(a) if \s z \ > sin(a), 



and instead of (28) the following symbol P ~ H\ will be used 



P = 2^(,-^S(<aM>)Y (30) 

Cmean \ C{y, Z) Ul J 

The function 

%, Z) C,a;)~5(^) (31) 
c(y,z) uj 

is a symbol (here we ignore that S is only C 2 instead of C°° , because a smoother symbol can 
be designed if the numerical computations require it). We define an operator B(y,z, D z ) by 

B(y, z, D z )u(y, z) = (2n)-^ [ [ (B(y, z, () + l.o.t.)e"< u(y, Q d(. (32) 



For later reference we also define a pseudodifferential operator B(y, z, D z , Dt) that acts on 
function of (z,t), with t a time coordinate 

B(y, z, D z ,D t )u(y, z, t) = (2vr)" d / / (B(y, z, (, w) + Lo.t.)e^-^%, C, -w) d( doj. (33) 



In these equations u defines a Fourier transform w.r.t. a subset of the variables as indicated. 
We will assume the lower order terms are such that B(y,z,D z ) and B(y, z, D z , Dt) are real 
and selfadjoint, see e.g. The differential equation assocated with (pOl) is 



(d y -iB(y,z,D z ))u(y,z)=0. (34) 

It is an evolution equation in y. 



Equation (34) is of a well known type, it is called a one-way wave equation. The unknown 
u can be considered as a function of (y,z,oo). After Fourier transform —uj h-> t, this results 
in a first order hyperbolic pseudodifferential equation results, see e.g. |21[ |3"T] for the theory 
of such equations. One-way wave equations are a tool for the analysis of hyperbolic systems 
of equations, as well as for their numerical computation, see [26] and the references therein. 
Because of their nature, they eiter describe waves propagating in the positive y direction, 
or in the negative y direction. So called 'turning waves', for which the y-component of the 
propagation direction changes sign, are not described by such an equation. Propagation 
under large angles with the y-axis is also a problem for numerical methods for one-way wave 
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equations. These usually become inaccuarate if the angle is too large, e.g. larger than 60 
degrees. 

By E(y, j/q) we will denote a solution operator, meaning if uq = uq(z) are initial conditions 
at yo, then 

(d y -iB)E(y,y )u = 0, E(y ,y )u = u. (35) 
We define the operator T, acting on a function of (y, z), by 

Tu(y,-) = E(y,0)u(y,-) (36) 

From the assumption that the operator B = B(y, z, D z ) is self-adjoint it follows that this 
operator is unitary, T*T = I. We find that 

(d y - iB)E(y, 0)u(y, •) = {{d y - iB)E(y, 0)) u(y, •) + E(y, 0)d y u 

= E(y,0)d y u, 

i.e. the following intertwining operator relation holds 

(d y - iB)T = Td y . (38) 

Let 4>[i(y), <pu{z) be frames, and suppose the frame ^(y) is adapted to the operator d y + 1 
in the sense discussed above, with preconditioning factors w^. We define a candidate for our 
frame by 

<?W(y> z ) = r< My)<M z )- ( 39 ) 

Since T is unitary the new frame is also a tight frame with frame bound 1. Indeed, 

Ei^>/)i 2 = Ei^/)i 2 = Ei^r*/)i 2 

MM A* (40) 

= (T*f,T*f) = 



We use the intertwining property and the unitarity property, to argue that the preconditioning 
weights should be w^. A matrix element satisfies 

(K,x,(dx-iB)^ u ) = (0 K 4>x,T*(d y -iB)T^4> u ) = (4,^)(0 A ,^), (41) 



hence, taking into account the factor — — in (28), Cl S ca " w u is a good choice for the precondi- 
tioning factors. 



Actual application of the frame transform associated with (39) requires the computation 
of E(y, 0)4> u (z). This amounts to solving the initial value problem for a first order hyperbolic 
pseudodifferential equation. There are various methods to solve one-way wave equations 
numerically, see [26j and references. Here we will use a different, semi-analytical approach, 
and use the WKB approximation. It is assumed that the frame property of the <j)^ tU , and the 
property that it yields is a good preconditioner are stable against the (small) perturbations 
due to the errors in the WKB approximation. 

For E(y, 0)(j) u (z) we therefore look for solutions of the form 

A{y,z)e i<t,{y ' z) . (42) 
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The "initial values" 4> u must be of the form A(z)e l ^ z \ Making use of our periodic setting, 
we choose a Fourier basis for <j) v (here v = k & Z) 

^ k = L~ l/2 e ix ^ k/L) . (43) 

An alternative would be a windowed Fourier frame 4>k, m = w {y ~ kvo)e irn< * m . 

Let's start with the equation for <j>. We write <j) = ujT(y, z), where T has the dimension of 
time. The initial conditions for T are T(0,z) = ^t-z. The WKB method gives the following 
equation for T 

dT c dT 

IT = (44) 
ay oj dz 



Without the regularization of (|29|) this becomes a well-known equation, 

0, (45) 



dT 1 (dT\ 2 



dy V c 2 \dz J 



the eikonal equation as an evolution equation in the y direction. Equation (44) can be solved 
e.g. by upwind finite differences. 

The amplitudes for the equation satisfy a transport equation along the characteristics |15j . 
Let Z(y, z) describe a set of characteristics. The fact that E is unitary provides an additional 
relation, namely that 



A(y, Z(y, z )) = A(0, z ) l—j . (46) 

Instead of computing the characteristics, it is convenient to compute Zo(y,z), the inverse of 
the map zq h-> Z(y, zq). This function satisfies itself a transport equation that can be solved 
along with T(y, z). 

The last step in describing the frame functions is the choice of (j)^{y)- We choose again 
the Fourier modes, \x = j, 4>j{y) = e t2lT ^^ 2Ll \ 

Having described the frame functions, we look at the resulting transform, which, we recall 
is defined by taking 



{(f>ti,u, f) = JJ <?W(^ z )f(v, z ) dz dy- (47) 
This amounts to two steps. First the map 

Fiu(y, C fc ) = / A(y, z)e^ T ^'^u{y, z) dz. (48) 



This is a Fourier integral operator. Then the map 

cj,k = T x ^j(F x u{yXk))- (49) 

We will call the frame transform a Lagrangian wave packet transform because of the tiling, 
which, as discussed below, is associated with a canonical relation that is by definition a 
Lagrangian manifold. 

We have already discussed that, if u is considered as a function of (y,z,uj), and after 
inverse Fourier transform — u i— > t, in short, in the time domain, equation ( |34[ ) becomes a 
first order hyperbolic pseudodifferential equation. In the time domain the operator T is a 
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Fourier integral operator. This is useful, because the theory of Fourier integral operators also 
includes a geometrical description of the mapping of energy in the phase space, which results 
(for large ui) in the description of the mapping of singularities. The mapping of singularities 
is according to the canonical relation |XQ^ |22| [331 15T] . The canonical relation for the solution 
operator to first order hyperbolic systems is described in terms of the bicharacteristics, i.e. 
the characteristics of the eikonal equation. We use the notation B(y,z,£) for the (principal) 
symbol of the pseudodifferential operator B. The bicharacteristics are described by the system 
of ODE's 



dz _ B 
dy C 
dC 
dy 



B 



Based on the unregularized operator (28), the explicit equations are 



dz 
dy 

d( 
dy 



uj 2 dc 
c 3 dz 



oo 



-1/2 



-1/2 



Denote by (Z(y, Zq, Co), &(y, zq, Co)) the solution operator that describes the flow with initial 
values (zo, Co) at y = 0. Let 7^ be the tiling associated with and T v be the tiling associated 
with (j) u . This results in the following description of the tiling 
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{(y,Z(y,z , Co); r] + B(y,Z,@),Q(y,z ,(o)) ! (^o,Co) G %, (y,7]) G %}. 



(50) 



5 Implementation aspects 

In this section we discuss some aspects of the implementation of the frame transform and the 
preconditioner outlined in the previous section. The two main steps to be performed for the 
frame transform are the localization in phase space in scale and in angle described in the first 



part of this section, and the transform described in (48) and (49). We will comment on the 
discretization of each of those steps, and then say a few words on the structure of a program 
implementing the preconditioner. 



The localization in scale in the Fourier domain was described in the text around (27). This 
is straightforward to discretize. We choose k\ = 0.2w / c ma x> k 2 = w/c max , k 3 = 1.0w/c min, 
Cmin- The algorithm appears to be insensitive to the detailed values. In connection 
with the subsampling step described below, it is advantageous to use smaller values of k± if 
possible. For the windowing in angle we use iV a = 8 angles, —225, —180, . . . , 90 degrees. Three 
subsequent angles are used to define a cutoff function in angle, using appropriately translated 
and dilated versions of the standard cutoff h. 

After windowing in the Fourier domain a subsampling step is done. After localization, the 
Fourier transform of the signal is supported in a small subdomain of the original domain. A 
rectangle is taken that contains this subdomain, and the Fourier grid within the rectangle is 
used for the further calculations. This results in a function defined on a substantially coarser 
grid in space. 
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The final step associated with the restriction in the phase space is the rotation to the 
new spatial coordinates (y, z), and a resorting of the data into periodic bands. The rotation 
amounts to interpolation to a rotated grid. This is done by shear rotation which is an exactly 
invertible operation. It changes the grid parameters, but that is not a problem. The shear 
rotation is combined with the resorting of the data in periodic bands. The details of the 
bands depends on the angle. The total number of bands will be denoted by N^. The bands 
are chosen overlapping, and a cutoff function in the y coordinate is applied. 

For the application of the Fourier integral operator, first a preparation step needs to be 
done computing T(y,z,p z ) and A(y,z,p z ). The travel time T(y,z,p z ) is computed using 
upwind finite differences. These are also used to compute Zo(y, z,p z ). A is then computed 
from (46). The storage for A(yXk,z)e lu)T ^ y '^ k ^' y,) is the largest storage requirement for the 
algorithm. For an d dimensional domain with discretization size N in each direction, and 
subsampling toward f3N in each direction, this is about 2iV a (/3A r ) 3 . Some numbers from the 
numerical experiments are given in the next section. 

Application of the Fourier integral operator in ( 48 ) is done simply by replacing the integral 
by a summation. We have experimented with a partial implementation of the scheme described 
in [3], but there were no significant gains, this was partially related to the large subsampling 
we did before computing (48). 

This leads to the following program structure. First all preparatory tasks are performed 
by a routine prepare_helmprec. The main steps are two call the two routines pre- 
pare_filter_band and prepare.lwpt, that prepare the filtering and banding step, respec- 
tively the Lagrangian wave packet transform. The results are stored and make it possible to 
apply the preconditioner rapidly each time it is called. The actual execution of the precon- 
ditioner is done by a routine execute_helmprec. A similar subdivision is made, a routine 
execute_filter_band executes the transform from input data to the filtered and banded 
data, and the adjoint of this operation. A routine executejlwpt computes the forward 
version or the adjoint of the transform described in (48) and (49). It is straightforward to 
apply the scaling factors u> M . 



6 Numerical results 

The purpose of our numerical experiments is to study the convergence of the method. The 
convergence depends on the problem parameters, the function c, the frequency cj, the damping 
parameter and the right hand side of the equation. Of this, the dependence on c and on u 
are the most interesting, and we will show some results for different values of c and oj. We 
have fixed the damping to a = 2ir. A smaller value leads to somewhat slower convergence, 
as one might expect. The right hand side will be chosen as a random array. We find that 
the number of steps to convergence depends very little on which random right hand side is 
chosen. 

The convergence speed is also influenced by the parameters in the algorithm. The dis- 
cretization was chosen to be 16-20 points per wavelength, where the minimum wavelength 
present in the domain was used. The discretization needs to be fine enough to have accurate 
dispersion relation, a finer discretization is not important for the convergence. We used 8 
angular directions, and two bands per angle. The number of bands can increase due to a 
combination of large propagation distances (i.e. width of the band in y direction) and strong 
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Figure 2: Some examples for the velocity c 



gradients, in which case the WKB method might no longer be valid. We did not observe 
this with the class of media chosen. After the filtering in scale and angle, the computations 
were done on a coarse grid of about two points per wavelength, which lead to a large re- 
duction in cost compared to the original grid. The number of iterations will be determined 
by requiring that the error is reduced with a factor of 10~ 5 in the LSQR method using the 
right-preconditioned form in ([6]). The convergence was observed to be linear, which gives an 
indication for the other error reduction factors. Other iterative methods were tried, such as 
BiCG and BiCGStab, but those were found to perform considerably worse. This is probably 
related to the distribution of the eigenvalues in the domain, eigenvalues were present near the 
positive and near the negative real axis. 

To study the dependence on c, we let it vary within a class of smooth functions, consisting 
of a sum of a few sines and cosine^] The domain will be the unit square. Some examples 
for c are given in Figure [2] A set of 100 computations were performed with different c and 
uj = 40-7T, corresponding to about 20 wavelengths in the domain. A histogram of the results 
is given in Figure [3j There is indeed some spread in the convergence speed. In general the 
media with higher velocity contrasts require more iterations, although we have not found a 
precise relation. 

Our main result concerns the dependence of the convergence on uj. In Figure [4] the number 
of iterations as a function of uj is given for a few choices of c. The main conclusion is that 
the number of iterations tends to become constant for the larger values of uj. For the smaller 
values of uj the convergence is slightly faster, by one or a few iterations. The results of Figure [4] 
were quite typical as c was varied, as can be seen from the comparison with Figure [3] 



7 Conclusion and discussion 

We presented a preconditioning method for the multi-dimensional Helmholtz equation. The 
preconditioning was based on a frame of functions designed to have a special phase space 
tiling, adapted to the Helmholtz operator. The tiles follow to a certain degree the level sets 

1 To be precise, this set is of the form 

c(x) = 1 + Re ^Ti.ifo +ib i )e 27ri{l ^ 2X1+1 '- 3X2)SS j (51) 

where aj,bj are random, uniformly distributed in [—1,1], and the triples jj take the values (0.12,1,0), 
(0.12,0,1), (0.084,1,1), (0.084,1,-1), (0.06,2,0), (0.06,0,2). 
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Iterations required: histogram 



15 20 25 

Niter 



Figure 3: Histogram describing the number of iterations for convergence for oj = 40ir and 100 
different choices of c. 



Iterations required vs. frequency for several examples 




100 200 300 400 500 
omega 



Figure 4: Number of iterations required vs. frequency for four choices of c within the class 
described 



of the absolute value of the symbol, they are of the general form 

|J {x} x E(x). (52) 

x&V 

Here is an rr-dependent subset of the Fourier coordinate space, V we chosen as a coor- 
dinate block. The tiles were curved because of the nontrivial dependence of on x. 

Our main goal was to establish whether this kind of phase space tiling could be useful 
for preconditioning the Helmholtz equation. This is clearly shown by the numerical experi- 
ments. Apart from the convergence, a second important property was that a large part of the 
computations could be done on very coarse grids. 

The generalization to 3-D would of course be of interest. A question is whether in 3- 



D the evaluation of the FIO in (48) can be speeded up by using the results of It is 
perhaps useful to point out some other directions for further development. One issue is the 
inclusion of boundary conditions. For applications like seismic imaging, the use of a simple 
absorbing boundary layers can be sufficient. The next step would be to include a boundary 
with Dirichlet or Neumann conditions at a planar or non-planar surface. Less smooth media 
are another interesting issue. We believe the results of this paper form a strong motivation 
for such further research. 
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